This invention relates generally to the field of seismic prospecting. More particularly, the invention relates to a method for constructing a model seismic image of a subsurface seismic reflector.
In the oil and gas industry, seismic prospecting techniques are commonly used to aid in the search for and evaluation of hydrocarbon deposits located in subterranean formations. In seismic prospecting, seismic energy sources are used to generate a seismic signal which propagates into the earth and is at least partially reflected by subsurface seismic reflectors. Such seismic reflectors typically are interfaces between subterranean formations having different elastic properties. The reflections are caused by differences in elastic properties, specifically wave velocity and rock density, which lead to differences in impedance at the interfaces. The reflections are recorded by seismic detectors at or near the surface of the earth, in an overlying body of water, or at known depths in boreholes. The resulting seismic data may be processed to yield information relating to the geologic structure and properties of the subterranean formations and potential hydrocarbon content.
The goal of all seismic data processing is to extract from the data as much information as possible regarding the subterranean formations. In order for the processed seismic data to fully represent geologic subsurface properties, the true amplitudes resulting from reflection of the input signal by the geologic target must be accurately represented. This requires that the amplitudes of the seismic data must be processed free from non-geologic seismic effects. Non-geologic amplitude effects include mechanisms that cause the measured seismic amplitudes to deviate from the amplitude caused by the reflection coefficient of the geologic target. These non-geologic amplitude effects can be related to acquisition of the data or to near surface effects. Examples of non-geologic amplitude effects that can be particularly troublesome are source and receiver variation, coherent noises, electrical noise or spikes, and overburden and transmission effects. If uncorrected, these effects can dominate the seismic image and obscure the true geologic picture.
A seismic wave source generates a wave that reflects from or xe2x80x9cilluminatesxe2x80x9d a portion of a reflector. The collection of sources that comprises an entire 3-D survey generally illuminates a large region of the reflector. Conventional prestack 3-D migration algorithms can produce precise images of the reflector only if illumination is relatively uniform. Lateral velocity variations within the earth and nonuniformly sampled 3-D prestack seismic data, however, generally cause reflectors to be illuminated nonuniformly. Nonuniform illumination is generally due to varying azimuths and source-receiver midpoint locations. Consequently, prestack 3-D migrated images are often contaminated with non-geologic artifacts called an xe2x80x9cacquisition footprintxe2x80x9d. These artifacts can interfere with the ultimate interpretation of seismic images and attribute maps. Understanding and removing the effects of the acquisition footprint has thus become important for seismic acquisition design, to processing, and interpretation.
Let {right arrow over (S)}=(S1, S2)T and {right arrow over (G)}=(G1, G2)T denote two-dimensional coordinate vectors of a seismic source, commonly called a shot, and a seismic receiver, typically a geophone, respectively. These vectors {right arrow over (S)} and {right arrow over (G)} are defined with respect to a global Cartesian coordinate system {right arrow over (x)}=(x, y, z) in the plane given by z=0. Schleicher et al., xe2x80x9c3-D True-Amplitude Finite-Offset Migrationxe2x80x9d, Geophysics, 58, 1112-1126, (1993) showed that for any specified measurement configuration of sources and receivers, the vector pair ({right arrow over (S)},{right arrow over (G)}) can be described by a single position vector {right arrow over ("xgr")}=("xgr"1, "xgr"2)T according to the following relations
{right arrow over (S)}({right arrow over ("xgr")})={right arrow over (S)}0+xcex93S({right arrow over ("xgr")}xe2x88x92{right arrow over ("xgr")}0)
and
{right arrow over (G)}({right arrow over ("xgr")})={right arrow over (G)}0+xcex93G({right arrow over ("xgr")}xe2x88x92{right arrow over ("xgr")}0).
Here {right arrow over (S)}0 and {right arrow over (G)}0 are coordinate vectors describing a fixed source-receiver pair ({right arrow over (S)}0, {right arrow over (G)}0) defined by position vector {right arrow over ("xgr")}={right arrow over ("xgr")}0. Configuration matrices xcex93S and xcex93G are 2xc3x972 constant matrices, depending only upon the measurement configuration. The configuration matrices are determined by             Γ      Sij        =                            ∂                      S            i                                    ∂                      ξ            j                              ⁢                        xe2x80x83                ⁢                  xe2x80x83                    ⁢      and                          Γ                  G          ij                    =                        ∂                      G            i                                    ∂                      ξ            j                                ,                  xe2x80x83            ⁢              xe2x80x83              ⁢                  for        ⁢                  xe2x80x83                ⁢        i            =      1        ,                  2        ⁢                  xe2x80x83                ⁢        and        ⁢                  xe2x80x83                ⁢        j            =      1        ,    2.  
Examples of sets of configuration matrices for particular measurement configurations are:
xcex93S=I and xcex93G=I for common offset configuration,
xcex93S=0 and xcex93G=I for common shot configuration, and
xcex93S=I and xcex93G=0 for common receiver configuration.
Here I is the 2xc3x972 identity matrix and 0 is the 2xc3x972 zero matrix.
As discussed by Schleicher et al., (1993), supra, common offset 3-D migration can be formulated as a weighted summation along the diffraction traveltime surface, xcfx84D, also called the Huygens surface, through the following Kirchhoff migration equation                               U          ⁡                      (                          x              →                        )                          =                                            ∑              j                        ⁢                          Δ              ⁢                              xe2x80x83                            ⁢                                                ξ                  →                                j                            ⁢                              w                ⁡                                  (                                                                                    ξ                        →                                            j                                        ,                                          x                      →                                                        )                                            ⁢                                                U                  .                                ⁡                                  (                                                                                    ξ                        →                                            j                                        ,                                          t                      +                                                                        τ                          D                                                ⁡                                                  (                                                                                                                    ξ                                →                                                            j                                                        ,                                                          x                              →                                                                                )                                                                                                      )                                                              ⁢                      ❘                          t              =              0                                .                                    (        1        )            
Here U is the seismic amplitude sum; {right arrow over (x)} is the image point; {right arrow over ("xgr")}j is the source-receiver midpoint for trace j; w is a weight chosen to preserve seismic amplitudes; and xcfx84D is the diffraction traveltime connecting source S, image point Q, and receiver G. The dot above U on the right hand side of equation (1) represents time differentiation. The term xcex94{right arrow over ("xgr")}j in migration equation (1) is the two-dimensional element of surface area representing trace j, indicating that this equation is a discretized integral formula. Neglect of variations in xcex94{right arrow over ("xgr")}j, caused by variations in source-receiver midpoint distribution, is in practice a significant cause of acquisition footprints. Migration implemented by equation (1) is referred to in the literature as xe2x80x9cTrue Amplitude Kirchhoff Migrationxe2x80x9d.
Consider an image volume consisting of N points in each of the x, y, and z directions. Then imaging with migration equation (1) is an O(N5) process, because the two-dimensional summations over {right arrow over ("xgr")}j must be performed at each point {right arrow over (x)} of the three-dimensional image space. O(N5) processes such as migration equation (1) must usually be implemented on supercomputers or on parallel networks of fast workstations.
Migration equation (1) only produces a good image under the assumptions that the source-receiver offset and azimuth are fixed in magnitude within a narrow range, and the data coverage is dense in {right arrow over ("xgr")}j-space. This last assumption means that there are no large data gaps, which would give large xcex94{right arrow over ("xgr")}j values. Failure of the data or the migration implementation to meet these assumptions generally results in an image that is contaminated with an acquisition footprint.
Schleicher et al., (1993), supra, show that synthetic seismic data may be written for one reflector and a collection of sources and receivers as the following data equation
U({right arrow over ("xgr")}j, t)=R({right arrow over ("xgr")}j)g(txe2x88x92xcfx84R({right arrow over ("xgr")}j))/xcex9({right arrow over ("xgr")}j).xe2x80x83xe2x80x83(2)
Here R is the reflection coefficient for the reflector, g(t) represents the seismic wavelet, xcfx84R is the reflection traveltime, and xcex9 is the geometrical spreading. If there were more than one reflector, an additional factor to account for the total loss in amplitude due to transmissions across the interfaces along the reflection raypath would be needed. Migration equation (1) becomes, upon insertion of synthetic data equation (2),                               U          ⁡                      (                          x              →                        )                          =                              ∑            j                    ⁢                      Δ            ⁢                          xe2x80x83                        ⁢                                          ξ                →                            j                        ⁢                          xe2x80x83                        ⁢                                                                                                      w                      ⁡                                              (                                                                                                            ξ                              →                                                        j                                                    ,                                                      x                            →                                                                          )                                                              ⁢                                          R                      ⁡                                              (                                                                              ξ                            →                                                    j                                                )                                                              ⁢                                                                  g                        .                                            ⁡                                              (                                                  t                          +                                                                                    τ                              D                                                        ⁡                                                          (                                                                                                                                    ξ                                    →                                                                    j                                                                ,                                                                  x                                  →                                                                                            )                                                                                -                                                                                    τ                              R                                                        ⁡                                                          (                                                                                                ξ                                  →                                                                j                                                            )                                                                                                      )                                                                              ⁢                                      ❘                                          t                      =                      0                                                                                        Λ                  ⁡                                      (                                                                  ξ                        →                                            j                                        )                                                              .                                                          (        3        )            
The target-oriented imaging equation (3) can be implemented as an O(N4) process because the image location {right arrow over (x)} may be restricted to only those points in image space that lie on the reflector surface.
In Kirchhoff migration in general, the set of {right arrow over ("xgr")} values summed over is referred to as the xe2x80x9cmigration aperturexe2x80x9d. Consider imaging equation (3) for some image point or reflection point {right arrow over (x)} on the reflector. Schleicher et al., xe2x80x9cMinimum Apertures and Fresnel Zones in Migration and Demigrationxe2x80x9d, Geophysics, 62, 183-194, (1997), showed that the optimum migration aperture contains only those traces that are reflected from the Fresnel zone area surrounding image point {right arrow over (x)}. This result only applies for migrating synthetic data that are free of diffracting targets, since otherwise a large aperture is needed to adequately focus diffracted energy.
The formulas derived by Schleicher et al. (1993) and (1997) apply only to an idealized situation in which offset and azimuth are strictly constant and source-receiver spacing is uniform and dense. In addition, their discussion of Fresnel zones is in relation to theoretical considerations about migration aperture. They did not consider the idea, introduced here for the present invention, of using weighted sums of Fresnel zone amplitudes to construct model seismic images which respect the exact field acquisition geometry in which source-receiver offset and azimuth vary over the survey.
The paper xe2x80x9cQuantifying Seismic Amplitude Distortions Below Saltxe2x80x9d by D. R. Muerdter, R. O. Lindsey, and D. W. Ratcliff: 1997 Offshore Technology conference paper OTC 8339, discusses the effect of irregular reflector illumination. The authors sort the results of raytrace modeling into reflection-point gathers and then prepare and analyze reflection point bin maps around salt structures. However, this method is only an approximation of seismic amplitudes.
Thus, a need exists for a method of removing non-geologic amplitude effects without adversely affecting the quality of the resulting data. There is a need for a method to reveal the acquisition footprint by modeling the migrated image more accurately, more economically, and faster than can be achieved by direct implementation of migration equation (1). There is a need for a method that creates maps of the acquisition footprint that are consistent with known wave theory of how seismic images are formed.
The present invention is a method for constructing a model seismic image of a subsurface seismic reflector. First, a set of source and receiver pairs is located and a subsurface velocity function is determined. Specular reflection points are determined on the subsurface seismic reflector for each of the source and receiver pairs. A Fresnel zone is determined on the subsurface seismic reflector for each of the specular reflection points, using the subsurface velocity function. One or more seismic wavelets are selected. A set of image points is defined containing the subsurface seismic reflector. A synthetic seismic amplitude is determined at each of the image points by summing the Fresnel zone synthetic seismic amplitudes for all of the Fresnel zones that contain the image point, using the seismic wavelets. Finally, the model seismic image of the subsurface seismic reflector is constructed, using the synthetic seismic amplitudes at the image points.